find.ladder <- function (x, ladder, ci.upp = 1.96, ci.low = 1.96, draw = TRUE,
                         dev = 50, warn = TRUE, init.thresh = 250, sep.index = 8,
                         method = "iter", avoid = 1000, who = "sample", attempt = 10,
                         cex.title = 0.8)
{
  MSE <- function(x, y) {
    X <- cbind(1, x)
    qr.X <- qr(X)
    b <- t(qr.Q(qr.X)) %*% y
    R <- qr.R(qr.X)
    beta <- as.vector(backsolve(R, b))
    fit <- X %*% beta
    res <- list(mse = sum((y - fit)^2), beta = beta[2])
    return(res)
  }
  MSE2 <- function(x, y) {
    X <- cbind(1, x)
    qr.X <- qr(X)
    b <- t(qr.Q(qr.X)) %*% y
    R <- qr.R(qr.X)
    beta <- as.vector(backsolve(R, b))
    fit <- X %*% beta
    mse <- sum((y - fit)^2)
    sst <- sum((y - mean(y))^2)
    r2 <- 1 - (mse/sst)
    res <- list(mse = mse, beta = beta, r2 = r2)
    return(res)
  }
  thresh = init.thresh
  roxy <- big.peaks.col(x[1:length(x)], thresh)
  nnn <- length(roxy$pos)
  while (nnn < length(ladder)) {
    cat("\nReducing threshold 2x to find ladder \n")
    thresh = thresh/2
    roxy <- big.peaks.col(x[1:length(x)], thresh)
    nnn <- length(roxy$pos)
  }
  whot <- length(roxy$pos) * 0.2
  what <- which(roxy$hei == max(roxy$hei))
  roxy <- separate(roxy, shift = sep.index, type = "pos")
  ii <- which(roxy$hei == max(roxy$hei)) + 1
  iii <- length(roxy$hei)
  roxy <- list(pos = roxy$pos[ii:iii], hei = roxy$hei[ii:iii])
  if (method == "iter") {
    step1 <- combn(roxy$pos[1:attempt], 3)
    step2 <- apply(step1/10, 2, MSE, y = ladder[1:3])
    mse <- unlist(lapply(step2, function(x) {
      x$mse
    }))
    covs <- apply(step1, 2, function(x, y) {
      cov(x, y)
    }, y = ladder[1:3])
    step2 <- mse * covs
    step3 <- step1[, which(step2 < sort(step2, decreasing = FALSE)[20])]
    step4 <- apply(step3, 2, function(x, y) {
      which(y %in% x)
    }, y = roxy$pos)
    caller <- function(roxy, www, ladder.call, x) {
      threshold <- length(x)
      posi <- numeric()
      fact2 <- length(ladder.call)
      expect <- roxy$pos[www]
      xxx <- ladder.call[c(1:3)]
      modx <- lm(expect ~ poly(xxx, degree = 1))
      expecto <- predict(modx, data.frame(xxx = ladder.call))
      ladder.call <- ladder.call[which(expecto < threshold * 0.85)]
      available <- length(roxy$pos) - length(ladder.call)
      ava2 <- length(ladder.call) - abs(available)
      if (available > 0) {
        if ((length(ladder.call) - 1) < 3) {
          tope <- length(ladder.call)
        } else {
          tope <- length(ladder.call) - 1
        }
      } else {
        tope <- ava2 - 2
      }
      expect <- rep(NA, tope + 1)
      for (i in 3:tope) {
        if (i == 3 & i != tope) {
          expect[1:3] <- roxy$pos[www]
          xxx <- ladder.call[c(1:3)]
          mod <- MSE2(xxx, expect[1:3])
          beta <- (mod)$beta
          expecto <- as.vector(beta[1] + matrix(ladder.call) %*% beta[-1])
          act <- roxy$pos[-which(roxy$pos %in% expect)]
          yoyo <- abs(expecto[i + 1] - act)
          good <- which(yoyo == min(yoyo))
          expect[i + 1] <- act[good]
          if (mod$r2 < 0.9) {
            i = tope
          }
        }
        if (i > 3 & i <= 5) {
          xx <- ladder.call[c(1:i)]
          mod <- MSE2(xx, expect[1:i])
          beta <- (mod)$beta
          expecto <- as.vector(beta[1] + matrix(ladder.call) %*% beta[-1])
          act <- roxy$pos[-which(roxy$pos %in% expect)]
          yoyo <- abs(expecto[i + 1] - act)
          good <- which(yoyo == min(yoyo, na.rm = TRUE))
          expect[i + 1] <- act[good]
          if (mod$r2 < 0.9) {
            i = tope
          }
        }
        if (i > 5) {
          xx <- cbind(ladder.call[c(1:i)], ladder.call[c(1:i)]^2,
                      ladder.call[c(1:i)]^3, ladder.call[c(1:i)]^4)
          mod <- MSE2(xx, expect[1:i])
          beta <- (mod)$beta
          if (length(which(is.na(beta))) > 0) {
            beta[which(is.na(beta))] <- 0
          }
          toto <- cbind(matrix(ladder.call), matrix(ladder.call)^2,
                        matrix(ladder.call)^3, matrix(ladder.call)^4)
          expecto <- cbind(rep(1, dim(toto)[1]), toto) %*% beta
          act <- roxy$pos[-which(roxy$pos %in% expect)]
          yoyo <- abs(expecto[i + 1] - act)
          good <- which(yoyo == min(yoyo))
          expect[i + 1] <- act[good]
          if (is.na(mod$r2)) {
            mod$r2 <- 0.1
          }
          if (mod$r2 < 0.9) {
            i = tope
          }
        }
        if (i == tope & i != 3) {
          if (i < 5) {
            expect[1:3] <- roxy$pos[www]
            xx <- ladder.call[c(1:i)]
          } else {
            xx <- cbind(ladder.call[c(1:i)], ladder.call[c(1:i)]^2,
                        ladder.call[c(1:i)]^3, ladder.call[c(1:i)]^4)
          }
          mod <- MSE2(xx, expect[1:i])
          beta <- (mod)$beta
          if (length(which(is.na(beta))) > 0) {
            beta[which(is.na(beta))] <- 0
          }
          if (i < 5) {
            toto <- cbind(matrix(ladder.call))
          } else {
            toto <- cbind(matrix(ladder.call), matrix(ladder.call)^2,
                          matrix(ladder.call)^3, matrix(ladder.call)^4)
          }
          expecto <- cbind(rep(1, dim(toto)[1]), toto) %*% beta
          act <- roxy$pos[-which(roxy$pos %in% expect)]
          yoyo <- abs(expecto[i + 1] - act)
          good <- which(yoyo == min(yoyo))
          expect[i + 1] <- act[good]
        }
        if (i == tope & i == 3) {
          expect[1:3] <- roxy$pos[www]
        }
      }
      posi <- expect
      tutu <- abs(length(x) - posi)
      posi <- posi[1:which(tutu == min(tutu, na.rm = TRUE))]
      heii <- roxy$hei[which(roxy$pos %in% posi)]
      fact3 <- length(posi)/fact2
      if (length((posi)) < 6) {
        fact <- summary(lm(ladder.call[1:length(posi)] ~
                             poly(posi, degree = length((posi)) - 1)))$r.squared * fact3
      } else {
        fact <- summary(lm(ladder.call[1:length(posi)] ~
                             poly(posi, degree = 5)))$r.squared * fact3
      }
      roxy2 <- list(pos = posi, hei = heii, wei = ladder.call[1:length(posi)],
                    corr = abs(cor(ladder.call[1:length(posi)], posi)),
                    error = fact)
      return(roxy2)
    }
    rt <- apply(data.frame(step4), 2, FUN = caller, roxy = roxy,
                ladder.call = ladder, x = x)
    corrs3 <- unlist(lapply(rt, function(x) {
      x$error
    }))
    roxy3 <- rt[[which(corrs3 == max(corrs3))]]
    if (draw == TRUE) {
      limi <- sort(roxy3$hei, decreasing = TRUE)
      plot(x, type = "l", xaxt = "n", ylim = c(0, (limi[3] + 1000)),
           cex.axis = 0.6, las = 2, xlim = c((min(roxy3$pos) - 100), (max(roxy3$pos) + 100)),
           col = transp("grey35", 0.7), ylab = "RFU", xlab = "",
           lwd = 2, main = attributes(x)$mycomm, cex.main = cex.title)
      axis(1, at = roxy3$pos, labels = roxy3$wei, cex.axis = 0.6)
      points(x = roxy3$pos, y = roxy3$hei, cex = 1.1, col = transp("black", 0.85))
      points(x = roxy3$pos, y = roxy3$hei, pch = 20, col = transp("red", 0.7))
      legend("topleft", legend = paste("Correlation:", round(roxy3$corr, digits = 4), sep = ""),
             bty = "n")
      legend("topright", legend = c("Peaks selected"),
             col = c("red"), bty = "n", pch = c(20), cex = 0.85)
    }
    roxy <- roxy3
  }
  if (method == "cicor") {
    roxy <- separate(roxy, shift = 40, type = "pos")
    nnn <- length(ladder)
    mm <- median(roxy$hei)
    vvv <- which(roxy$hei == max(roxy$hei))
    se <- sd(roxy$hei[-vvv])/sqrt(length(roxy$hei[-vvv]))
    reduced <- roxy$pos[which(roxy$hei < mm + (1 * se) &
                              roxy$pos > mm - (1 * se) &
                              roxy$hei > init.thresh)]
    reduced2 <- roxy$pos[which(roxy$pos >= min(reduced))]
    if (length(reduced2) >= nnn & length(reduced2) < nnn + 8) {
      all.combs <- combn(reduced2, m = length(ladder))
      cors <- apply(all.combs, 2, function(x, y) {
        cor(x, y)
      }, y = ladder)
      found <- all.combs[, which(cors == max(cors))]
      roxy <- list(pos = found, hei = x[found], wei = ladder)
      nono <- length(x) - max(roxy$pos)
      if (draw == TRUE) {
        plot(x, type = "l", xaxt = "n", ylab = "Intensity",
             col = transp("gray29", 0.6), main = attributes(x)$mycomm,
             cex.main = cex.title)
        axis(1, at = roxy$pos, labels = roxy$wei, cex.axis = 0.6)
        abline(v = found, col = "red", lty = 3)
      }
    } else {
      roxy <- list(pos = seq(1, nnn) + rnorm(nnn, 0, 1),
                   hei = seq(1, nnn) + rnorm(nnn, 0, 1),
                   wei = ladder)
      cat(file=stderr(), "The ladder has someting wrong! Maybe the experiment not be very well\n")
      #plot(x, type = "l", col = transp("black", 0.5),
      #     main = attributes(x)$mycomm, cex.main = cex.title)
    }
  }
  if (method == "red") {
    nnn <- length(roxy$pos) - length(ladder)
    if (nnn > 15) {
      yoy <- length(roxy$pos) - length(ladder)
      yoy2 <- rev(seq(1, yoy, by = 1))
      corres <- numeric()
      for (i in 1:yoy) {
        corres[i] <- cor(roxy$pos[i:(length(roxy$pos) - yoy2[i])], ladder)
      }
      vv <- which(corres == max(corres))[1]
      roxy <- list(pos = roxy$pos[vv:((vv - 1) + length(ladder))],
                   hei = roxy$hei[vv:((vv - 1) + length(ladder))],
                   wei = ladder)
      if (draw == TRUE) {
        plot(x, type = "l", xaxt = "n", ylim = c(0, (max(roxy$hei) + 1000)),
             cex.axis = 0.6, las = 2, xlim = c((min(roxy$pos) - 100), (max(roxy$pos) + 100)),
             col = transp("grey35", 0.7), ylab = "RFU",
             xlab = "", lwd = 2, main = attributes(x)$mycomm, cex.main = cex.title)
        axis(1, at = roxy$pos, labels = roxy$wei, cex.axis = 0.6)
        points(x = roxy$pos, y = roxy$hei, cex = 1.1, col = transp("black", 0.85))
        points(x = roxy$pos, y = roxy$hei, pch = 20, col = transp("red", 0.7))
        legend("topleft", legend = paste("Correlation found:", round(max(cors), digits = 4), sep = ""),
               bty = "n")
        legend("topright", legend = c("Peaks selected"),
               col = c("red"), bty = "n", pch = c(20), cex = 0.85)
      }
    } else {
      nnn <- length(ladder)
      le <- length(roxy$pos[1]:(roxy$pos[length(roxy$pos)]))
      tra <- (le * ladder)/max(ladder)
      tries <- length(seq(tra[length(tra)], roxy$pos[length(roxy$pos)], 1)) + 100
      abso <- function(test1, pos) {
        test1 <- matrix(test1, nrow = 1)
        res <- apply(test1, 2, function(x, y) {
          xx1 <- abs(as.vector(x) - y)
          z <- y[which(xx1 == min(xx1))[1]]
          return(z)
        }, y = pos)
        return(res)
      }
      all.tests <- apply(matrix(1:tries, nrow = 1), 2,
                         function(q1, q2) {
                           q1 + q2
                         }, q2 = tra)
      all.abso <- apply(all.tests, 2, abso, pos = roxy$pos)
      all.cors <- apply(all.abso, 2, function(x, y) {
        cor(x, y)
      }, y = ladder)
      step1 <- (all.tests - all.abso)^2
      all.ss2 <- apply(step1, 2, sum)
      sss2 <- abs((all.ss2 - mean(all.ss2))/sd(all.ss2))
      response <- all.cors/sss2
      vv4 <- which(response >= min(sort(response, decreasing = TRUE)[1:5]))
      reduced <- sort(unique(as.vector(all.tests[, vv4])))
      reduced2 <- roxy$pos[which(roxy$pos >= min(reduced))]
      if (length(reduced2) >= nnn & length(reduced2) < nnn + 8) {
        all.combs <- combn(reduced2, m = length(ladder))
        cors <- apply(all.combs, 2, function(x, y) {
          cor(x, y)
        }, y = ladder)
        found <- all.combs[, which(cors == max(cors))]
        roxy <- list(pos = found, hei = x[found], wei = ladder)
        if (draw == TRUE) {
          plot(x, type = "l", xaxt = "n", ylim = c(0, (max(roxy$hei) + 1000)),
               cex.axis = 0.6,
               las = 2, xlim = c((min(roxy$pos) - 100), (max(roxy$pos) + 100)),
               col = transp("grey35", 0.7), ylab = "RFU", xlab = "", lwd = 2,
               main = attributes(x)$mycomm, cex.main = cex.title)
          axis(1, at = roxy$pos, labels = roxy$wei, cex.axis = 0.6)
          points(x = roxy$pos, y = roxy$hei, cex = 1.1,
                 col = transp("black", 0.85))
          points(x = roxy$pos, y = roxy$hei, pch = 20,
                 col = transp("red", 0.7))
          legend("topleft", legend = paste("Correlation found:", round(max(cors), digits = 4), sep = ""),
                 bty = "n")
          legend("topright", legend = c("Peaks selected"),
                 col = c("red"), bty = "n", pch = c(20), cex = 0.85)
        }
      } else {
        roxy <- list(pos = seq(1, nnn) + rnorm(nnn, 0, 1),
                     hei = seq(1, nnn) + rnorm(nnn, 0, 1),
                     wei = ladder)
        print("Friend I don't think you have ladder in this sample or you used the wrong ladder, look the plot")
        plot(x, type = "l", main = attributes(x)$mycomm,
             cex.main = cex.title)
      }
    }
  }
  if (method == "cor") {
    if (length(roxy$pos) > (length(ladder) + 10)) {
      print(paste("WOOW too many peaks in this", who, "!! low thresholds throw too many noisy peaks, consider increasing the initial thresold for your ladder, the number of possible combinations is too high to be computed, we will have to do 15,000 samples and get the most likely sizing, you better double check this sample"))
      nono <- which(roxy$pos < avoid)
      roxy <- list(pos = roxy$pos[-nono], hei = roxy$hei[-nono])
      roxy <- separate(roxy, shift = sep.index, type = "pos")
      pos.mod <- matrix(0, ncol = 15000, nrow = length(ladder))
      for (k in 1:15000) {
        pos.mod[, k] <- sort(sample(roxy$pos, size = length(ladder), replace = FALSE),
                             decreasing = FALSE)
      }
      dd <- apply(pos.mod, 2, function(x5, ladder) {
        cor(x5, ladder)
      }, ladder)
      v <- which(dd == max(dd))[1]
      v2 <- which(roxy$pos %in% pos.mod[, v])
      roxy <- list(pos = pos.mod[, v],
                   hei = roxy$hei[v2],
                   wei = ladder)
    } else {
      mi = length(ladder)
      if (mi > length(roxy$pos)) {
        print(paste("ERROR!! using the initial threshold you specified we did not find enough peaks for",
                    who, ", please try reducing the 'init.thresh' or check the plot for this plant using 'detect.ladder' function, maybe this sample did not have ladder :)"))
        roxy <- list(pos = seq(1, length(ladder)) + rnorm(length(ladder), 0, 1),
                     hei = seq(1, length(ladder)) + rnorm(length(ladder), 0, 1),
                     wei = ladder)
      } else {
        pos.mod <- combn(roxy$pos, m = mi)
        dd <- apply(pos.mod, 2, function(x5, ladder) {
          cor(x5, ladder)
        }, ladder)
        v <- which(dd == max(dd))
        v2 <- which(roxy$pos %in% pos.mod[, v])
        roxy <- list(pos = pos.mod[, v],
                     hei = roxy$hei[v2],
                     wei = ladder)
      }
    }
  }
  if (method == "ci") {
    z <- which(roxy$hei == max(roxy$hei))
    mm <- median(roxy$hei[-z])
    se2.low <- (sd(roxy$hei[-z])/sqrt(length(roxy$hei[-z]))) * ci.low
    se2.upp <- (sd(roxy$hei[-z])/sqrt(length(roxy$hei[-z]))) * ci.upp
    v <- which((roxy$hei > (mm - se2.low)) & (roxy$hei < (mm + se2.upp)))
    roxy <- list(pos = roxy$pos[v], hei = roxy$hei[v])
    vv <- which(diff(roxy$pos) < dev)
    vv2 <- vv + 1
    while (length(vv) > 0) {
      keep <- numeric()
      for (h in 1:length(vv)) {
        a1 <- vv[h]
        a2 <- vv2[h]
        a3 <- c(roxy$hei[a1], roxy$hei[a2])
        a4 <- c(a1, a2)
        keep[h] <- (a4[which(a3 == max(a3))])[1]
      }
      keep <- unique(keep)
      "%!in%" <- function(x, y) !(x %in% y)
      keep2 <- unique(c(vv, vv2)[which(c(vv, vv2) %!in% keep)])
      roxy <- list(pos = roxy$pos[-keep2], hei = roxy$hei[-keep2])
      vv <- which(diff(roxy$pos) < dev)
      vv2 <- vv + 1
    }
    s1 <- length(roxy$pos) - (length(ladder) - 1)
    s2 <- length(roxy$pos)
    if (s1 <= 0) {
      if (warn == TRUE) {
        print("Are you sure this is a ladder channel, we did not find a clear pattern, stop a minute to check the plot")
      }
      if (draw != F) {
        plot(x, ylim = c(0, mm + se2.upp), type = "l",
             main = attributes(x)$mycomm, cex.main = cex.title)
      }
      roxy <- list(pos = seq(1, length(ladder)) + rnorm(length(ladder), 0, 1),
                   hei = seq(1, length(ladder)) + rnorm(length(ladder),
                                                                                                    0, 1), wei = ladder)
    } else {
      roxy <- list(pos = roxy$pos[s1:s2], hei = roxy$hei[s1:s2],
                   wei = ladder)
    }
  }
  if (method == "cor" | method == "ci") {
    if (draw == TRUE) {
      xx <- roxy$wei
      yy <- roxy$pos
      mod <- lm(yy ~ xx)
      xlabels <- as.vector(predict(mod, newdata = data.frame(xx = seq(0, max(ladder), by = 25))))
      plot(x, type = "l", main = attributes(x)$mycomm,
           cex.main = cex.title, xlab = "", ylab = "Intensity",
           xaxt = "n", col = transp("grey39", 0.6),
           xlim = c(0, (max(roxy$pos) + 100)),
           ylim = c(-100, max(roxy$hei)))
      if (method == "cor") {
        legend("topleft", legend = paste("Correlation=", round(max(dd)), sep = ""),
               bty = "n")
      }
      axis(1, at = xlabels, labels = seq(0, max(ladder), by = 25), cex.axis = 0.9)
      abline(v = roxy[[1]], col = "red", lty = 3)
      if (method == "ci") {
        abline(h = mm, col = "blue", lty = 3)
        abline(h = (mm + se2.upp), col = "blue", lty = 3)
        abline(h = (mm - se2.low), col = "blue", lty = 3)
        legend("topright", legend = c("90% CI", "Peaks selected"),
               col = c("blue", "red"), bty = "n", lty = c(3, 3), cex = 1)
      } else {
        legend("topright", legend = c("Peaks selected"),
               col = c("red"), bty = "n", lty = c(3), cex = 1)
      }
      text(x = roxy[[1]], rep(-200, length(ladder)), labels = ladder,
           cex = 0.6)
    }
  }
  return(roxy)
}
